home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / pi.cal < prev    next >
Text File  |  1995-07-17  |  1KB  |  55 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Calculate pi within the specified epsilon using the quartic convergence
  7.  * iteration.
  8.  */
  9.  
  10. define qpi(epsilon)
  11. {
  12.     local niter, yn, ym, tm, an, am, t, tn, sqrt2, epsilon2, count, digits;
  13.     local bits, bits2;
  14.  
  15.     if (isnull(epsilon))
  16.         epsilon = epsilon();
  17.     digits = digits(1/epsilon);
  18.     if    (digits <=  8) { niter = 1; epsilon =    1e-8; }
  19.     else if (digits <= 40) { niter = 2; epsilon =  1e-40; }
  20.     else if (digits <= 170) { niter = 3; epsilon = 1e-170; }
  21.     else if (digits <= 693) { niter = 4; epsilon = 1e-693; }
  22.     else {
  23.         niter = 4;
  24.         t = 693;
  25.         while (t < digits) {
  26.             ++niter;
  27.             t *= 4;
  28.         }
  29.     }
  30.     epsilon2 = epsilon/(digits/10 + 1);
  31.     digits = digits(1/epsilon2);
  32.     sqrt2 = sqrt(2, epsilon2);
  33.     bits = abs(ilog2(epsilon)) + 1;
  34.     bits2 = abs(ilog2(epsilon2)) + 1;
  35.     yn = sqrt2 - 1;
  36.     an = 6 - 4 * sqrt2;
  37.     tn = 2;
  38.     for (count = 0; count < niter; count++) {
  39.         ym = yn;
  40.         am = an;
  41.         tn *= 4;
  42.         t = sqrt(sqrt(1-ym^4, epsilon2), epsilon2);
  43.         yn = (1-t)/(1+t);
  44.         an = (1+yn)^4*am-tn*yn*(1+yn+yn^2);
  45.         yn = bround(yn, bits2);
  46.         an = bround(an, bits2);
  47.     }
  48.     return (bround(1/an, bits));
  49. }
  50.  
  51. global lib_debug;
  52. if (lib_debug >= 0) {
  53.     print "qpi(epsilon) defined";
  54. }
  55.